program drop _all
program define TABLE9prog, rclass
version 14
syntax, beta0(real) betax(real) betay(real) truLRP(real) numN(integer) numT(integer) ///
    rhox(real)
drop _all


// Create the data
set obs `numN'
gen id = _n
// "ai" is the part of the error term that doesn't change over time for a given
// unit.
gen ai = rnormal()
expand `numT'
bysort id: gen t=_n
xtset id t
gen x = 0
replace x = `rhox'*L.x + rnormal() if t > 1
// "uit" is the part of the error term that varies across time
gen uit = rnormal()
// The total error term is the sum of "a" and "e"
gen error = ai + uit
gen y = (`beta0')/(1-`betay')
replace y = `beta0' + `betax'*x + `betay'*L.y + error if t > 1
gen Ly = y
replace Ly = L.y if t > 1


////////////////////////////////////////////////////
// This section calculates DFE
////////////////////////////////////////////////////

// Wald
xtpmg d.y if t > 15, lr(l.y x) ec(ecdfe) dfe
test [ecdfe]x = `truLRP'
return scalar pDFEwald = r(p)

// Fieller
regress y x Ly i.id if t > 15
mat cov = e(V)
matrix VV = cov[1..2, 1..2]
scalar LRP0 = `truLRP'
scalar ztest = (_b[x]- LRP0*(1-_b[Ly]))/sqrt(LRP0^2*VV[2,2]+2*LRP0*VV[1,2]+VV[1,1])
return scalar pDFEfieller = (1 - normal(abs(ztest)))*2


////////////////////////////////////////////////////
// This section calculates DGMM
////////////////////////////////////////////////////

// Wald
xtabond y x if t > 15
testnl _b[x]/(1-_b[L.y])=`truLRP'
return scalar pDGMMwald = r(p)

// Fieller
mat cov = e(V)
matrix VV = cov[1..2, 1..2]
scalar LRP0 = `truLRP'
scalar ztest = (_b[x]- LRP0*(1-_b[L.y]))/sqrt(LRP0^2*VV[1,1]+2*LRP0*VV[1,2]+VV[2,2])
return scalar pDGMMfieller = (1 - normal(abs(ztest)))*2




////////////////////////////////////////////////////
// This section calculates SGMM
////////////////////////////////////////////////////

// Wald
xtdpdsys  y x if t > 15
testnl _b[x]/(1-_b[L.y])=`truLRP'
return scalar pSGMMwald = r(p)

// Fieller
mat cov = e(V)
matrix VV = cov[1..2, 1..2]
scalar LRP0 = `truLRP'
scalar ztest = (_b[x]- LRP0*(1-_b[L.y]))/sqrt(LRP0^2*VV[1,1]+2*LRP0*VV[1,2]+VV[2,2])
return scalar pSGMMfieller = (1 - normal(abs(ztest)))*2

end
